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1 Introduction 

Time-dependent density-functional theory (TDDFT) extends the basic ideas 
of ground-state density-functional theory (DFT) to the treatment of excita- 
tions or more general time-dependent phenomena. TDDFT can be viewed 
an alternative formulation of time-dependent quantum mechanics but, in 
contrast to the normal approach that relies on wave-functions and on the 
many-body Schrodinger equation, its basic variable is the one-body electron 
density, n(r,t). The advantages are clear: The many-body wave- function, a 
function in a 3iV-dimcnsional space (where N is the number of electrons in 
the system), is a very complex mathematical object, while the density is a 
simple function that depends solely on 3 variables, x, y and z. The standard 
way to obtain n(r, t) is with the help of a fictitious system of non-interacting 
electrons, the Kohn-Sham system. The final equations are simple to tackle nu- 
merically, and are routinely solved for systems with a large number of atoms. 
These electrons feel an effective potential, the time-dependent Kohn-Sham 
potential. The exact form of this potential is unknown, and has therefore to 
be approximated. 

The scheme is perfectly general, and can be applied to essentially any 
time-dependent situation. Two regimes can however be observed: If the time- 
dependent potential is weak, it is sufficient to resort to linear-response theory 
to study the system. In this way it is possible to calculate, e.g. optical absorp- 
tion spectra. It turns out that, even with the simplest approximation to the 
Kohn-Sham potential, spectra calculated within this framework are in very 
good agreement with experimental results. On the other hand, if the time- 
dependent potential is strong, a full solution of the Kohn-Sham equations is 
required. A canonical example of this regime is the treatment of atoms or 
molecules in strong laser fields. In this case, TDDFT is able to describe non- 
linear phenomena like high-harmonic generation, or multi-photon ionization. 

Our purpose, while writing this chapter, is to provide a pedagogical in- 
troduction to TDDFT 1 . With that in mind, we present, in section 2, a 
quite detailed proof of the Runge-Gross theorem[5], i.e. the time-dependent 
generalization of the Hohenberg-Kohn theorem [6], and the corresponding 

1 The reader interested in a more technical discussion is therefore invited to read 
Ref. [1-4], where also very complete and updated lists of references can be found. 
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Kohn-Sham construction [7]. These constitute the mathematical foundations 
of TDDFT. Several approximate exchange-correlation (xc) functionals arc 
then reviewed. In section 3 we are concerned with linear-response theory, 
and with its main ingredient, the xc kernel. The calculation of excitation en- 
ergies is treated in the following section. After giving a brief overlook of the 
competing density-functional methods to calculate excitations, we present 
some results obtained from the full solution of the Kohn-Sham scheme, and 
from linear-response theory. Section 5 is devoted to the problem of atoms 
and molecules in strong laser fields. Both high-harmonic generation and ion- 
ization are discussed. Finally, the last section is reserved to some concluding 
remarks. 

For simplicity, we will write all formulae for spin-saturated systems. Ob- 
viously, spin can be easily included in all expressions when necessary. Hartree 
atomic units (e = K = m = 1) will be used throughout this chapter. 

2 Time-dependent DFT 
2.1 Preliminaries 

A system of N electrons with coordinates r = (r\ ■ ■ ■ rjv) is known to obey 
the time-dependent Schrodinger equation 



This equation expresses one of the most fundamental postulates of quantum 
mechanics, and is one of the most remarkable discoveries of physics during the 
XXth century. The absolute square of the electronic wave- function, \&{r_, t)\ 2 , 
is interpreted as the probability of finding the electrons in positions r. The 
Hamiltonian can be written in the form 



i^-Mr,t) = H{r,tp(r,t) , 



(1) 



f(r) + W(r) + V C xt(L,t). 



(2) 



The first term is the kinetic energy of the electrons 




(3) 



while W accounts for the Coulomb repulsion between the electrons 





Furthermore, the electrons are under the influence of a generic, time-dependent 
potential, V cxt (r, t). The Hamiltonian (2) is completely general and describes 
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a wealth of physical and chemical situations, including atoms, molecules, and 
solids, in arbitrary time-dependent electric or magnetic fields, scattering ex- 
periments, etc. In most of the situations dealt with in this article we will be 
concerned with the interaction between a laser and matter. In that case, we 
can write the time-dependent potential as the sum of the nuclear potential 
and a laser field, Vtd = U cn + Vi asel . The term U CD accounts for the Coulomb 
attraction between the electrons and the nuclei, 

N n N 

V— 1 %— 1 1 v ' 1 

where Z v and R v denote the charge and position of the nucleus v, and N n 
stands for the total number of nuclei in the system. Note that by allowing the 
R v to depend on time we can treat situations where the nuclei move along a 
classical path. This may be useful when studying, e.g., scattering experiments, 
chemical reactions, etc. The laser field, Viaser, reads, in the length gauge, 

N 

Viaser (r, t) = E f(t) sin(wi) ^2 r i ■ a > ( 6 ) 

1=1 

where a, oj and E are respectively the polarization, the frequency and the 
amplitude of the laser. The function f(t) is a envelope that shapes the laser 
pulse during time. Note that in writing Eq. (6) we use two approximations: 
i) We treat the laser field classically, i.e., we do not quantize the photon field. 
This is a well justified procedure when the density of photons is large and the 
individual (quantum) nature of the photons can be disregarded. In all cases 
presented in this article this will be the case, ii) Expression (6) is written 
within the dipole approximation. The dipole approximation holds whenever 
(a) The wave-length of the light (A = 2itc/uj, where c is the velocity of light 
in vacuum) is much larger than the size of the system. This is certainly 
true for all atoms and most molecules we are interested in. However, one 
has to be careful when dealing with very large molecules (e.g. proteins) or 
solids, (b) The path that the particle travels in one period of the laser field 
is small compared to the wavelength. This implies that the average velocity 
of the electrons, v, fulfills vT < A ^ i; < A/T = c, where T stands for 
the period of the laser. In these circumstances we can treat the laser field 
as a purely electric field and completely neglect its magnetic component. 
This approximation holds if the intensity of the laser is not strong enough 
to accelerate the electrons to relativistic velocities, (c) The total duration of 
the laser pulse should be short enough so that the molecule does not leave 
the focus of the laser during the time the interaction lasts. 

Although the many-body Schrodinger equation, Eq. (1), achieves, to our 
best knowledge, a remarkably good description of nature, it poses a tanta- 
lizing problem to scientists. Its exact (in fact, numerical) solution has been 
achieved so far for a disappointingly small number of particles. In fact, even 
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the calculation of a "simple" two electron system (the Helium atom) in a laser 
field takes several months in a modern computer [8] (see also the work on the 
H^[9] molecule and the H^ + molecule [10]). The effort to solve Eq. (1) grows 
exponentially with the number of particles, therefore rapid developments re- 
garding the exact solution of the Schrodinger equation are not expected. 

In these circumstances, the natural approach of the theorist is to trans- 
form and approximate the basic equations to a manageable level that still 
retains the qualitative and (hopefully) quantitative information about the 
system. Several techniques have been developed throughout the years in the 
quantum chemistry and physics world. One such technique is TDDFT. Its 
goal, like always in density- functional theories, is to replace the solution of 
the complicated many-body Schrodinger equation by the solution of the much 
simpler one-body Kohn-Sham equations, thereby relieving the computational 
burden. 

The first step of any DFT is the proof of a Hohcnberg-Kohn type theo- 
rem^]. In its traditional form, this theorem demonstrates that there exists 
a one-to-one correspondence between the external potential and the (one- 
body) density. With the external potential it is always possible (in princi- 
ple) to solve the many-body Schrodinger equation to obtain the many-body 
wave-function. From the wave-function we can trivially obtain the density. 
The second implication, i.e. that the knowledge of the density is sufficient 
to obtain the external potential, is much harder to prove. In their seminal 
paper, Hohenberg and Kohn used the variational principle to obtain a proof 
by reductio ad absurdum. Unfortunately, their method cannot be easily gen- 
eralized to arbitrary DFTs. The Hohcnberg-Kohn theorem is a very strong 
statement: From the density, a simple property of the quantum mechani- 
cal system, it is possible to obtain the external potential and therefore the 
many-body wave-function. The wave-function, by its turn, determines every 
observable of the system. This implies that every observable can be written 
as a functional of the density. 

Unfortunately, it is very hard to obtain the density of an interacting sys- 
tem. To circumvent this problem, Kohn and Sham introduced an auxiliary 
system of non-interacting particles [7] . The dynamics of these particles are 
governed by a potential chosen such that the density of the Kohn-Sham 
system equals the density of the interacting system. This potential is local 
(multiplicative) in real space, but it has a highly non-local functional de- 
pendence on the density. In non-mathematical terms this means that the 
potential at the point r can depend on the density of all other points (e.g. 
through gradients, or through integral operators like the Hartrec potential). 
As we are now dealing with non-interacting particles, the Kohn-Sham equa- 
tions are quite simple to solve numerically However, the complexities of the 
many-body system are still present in the so-called exchange-correlation (xc) 
functional that needs to be approximated in any application of the theory. 



Time-dependent Density Functional Theory 



5 



2.2 The Runge-Gross theorem 

In this section, we will present a detailed proof of the Runge-Gross theorem[5] , 
the time-dependent extension of the ordinary Hohenberg-Kohn theorem [6]. 
There are several "technical" differences between a time-dependent and a 
static quantum-mechanical problem that one should keep in mind while try- 
ing to prove the Runge-Gross theorem. In static quantum mechanics, the 
ground-state of the system can be determined through the minimization of 
the total energy functional 



In time-dependent systems, there is no variational principle on the basis of 
the total energy for it is not a conserved quantity. There exists, however, a 
quantity analogous to the energy, the quantum mechanical action 



where <P(t) is a iVe-body function defined in some convenient space. From 
expression (8) it is easy to obtain two important properties of the action: 
i) Equating the functional derivative of Eq. (8) in terms of $*{t) to zero, we 
arrive at the time-dependent Schrodingcr equation. We can therefore solve the 
time-dependent problem by calculating the stationary point of the functional 
A[<P]. The function \P(t) that makes the functional stationary will be the 
solution of the time-dependent many-body Schrodingcr equation. Note that 
there is no "minimum principle" , but only a "stationary principle" . ii) The 
action is always zero at the solution point, i.e. A[&] = 0. These two properties 
make the quantum-mechanical action a much less useful quantity than its 
static counterpart, the total energy. 

Another important point, often overlooked in the literature, is that a 
time-dependent problem in quantum mechanics is mathematically defined as 
an initial value problem. This stems from the fact that the time-dependent 
Schrodingcr equation is a first-order differential equation in the time coordi- 
nate. The wave-function (or the density) thus depends on the initial state, 
which implies that the Runge-Gross theorem can only hold for a fixed ini- 
tial state (and that the xc potential depends on that state). In contrast, the 
static Schrodinger equation is a second order differential equation in the space 
coordinates, and is the typical example of a boundary value problem. 

From the above considerations the reader could conjecture that the proof 
of the Runge-Gross theorem is more involved than the proof of the ordinary 
Hohenberg-Kohn theorem. This is indeed the case. What we have to demon- 
strate is that if two potentials, v(r,t) and w'(r,t), differ by more thanapurely 
time dependent function c(i) 2 , they cannot produce the same time-dependent 

2 If the two potentials differ solely by a time-dependent function, they will produce 
wave-functions which are equal up to a purely time-dependent phase. This phase 



E[<P] = (<P\ H |<P) . 



(7) 
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density, n(r,t), i.e. 

v{r, t) ± v'{r, t) + c(t) => p(r, t) ^ p'(r, t) . 



(9) 



This statement immediately implies the one-to-one correspondence between 
the potential and the density. In the following we will utilize primes to dis- 
tinguish the quantities of the systems with external potentials v and v' . Due 
to technical reasons that will become evident during the course of the proof, 
we will have to restrict ourselves to external potentials that are Taylor ex- 
pandable with respect to the time coordinate around the initial time t 



v(r,t) = J2c k (r)(t-t ) k , 



k=0 



with the expansion coefficients 



1 d k 



t=t 



We furthermore define the function 

Qk 



u k {r)= Q^[v(r,t)-v'(r,t)] 



(10) 



(11) 



(12) 



t=t 



Clearly, if the two potentials are different by more than a purely time- 
dependent function, at least one of the expansion coefficients in their Taylor 
expansion around to will differ by more than a constant 



3fe>o : u k (r) ^ const. 



(13) 



In the first step of our proof we will demonstrate that if v ^ v' + c(t) then 
the current densities, j and j' , generated by v and v' are also different. The 
current density j can be written as the expectation value of the current 
density operator: 

j(r,t) = <*(*)! J'W !*(*)) . (14) 



where the operator j is written 

j(r) = ~ { [V^t(r)] ^(r) - ^(r) [v^(r)] } . 



(15) 



We now use the quantum-mechanical equation of motion, which is valid for 
any operator, 0(t), 



if t (*(t)\d(t)\*(t)) = (*(t)\i^6(t)+[d(t),H(t)] \nt)) 



(16) 



will, of course, cancel while calculating the density (or any other observable, in 
fact). 
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to write the equation of motion for the current density in the primed and 
unprimed systems 



i-j(r,t) =<*(*)! [j(r),H(t)\ |#(t)) 
i±j'(r,t) = (*'(t)\[j(r),H'(t)] \V'(t)) 



(17) 
(18) 



As we start from a fixed initial many-body state, at to the wave-functions, 
the densities, and the current densities have to be equal in the primed and 
unprimed systems 



|ff(to)) = \* '(to)) = Wo) 
n(r,t ) = n'(r,t ) = n (r) 
j{r,t ) =j'(r,t ) =j (r). 



(19) 
(20) 
(21) 



If we now take the difference between the equations of motion (17) and (18) 
we obtain, when t = t , 

i| [j(r,t) - j'(r,t)] t=tQ = (%\ [j(r),H(to) H'(t )} \%) 

= (%\ j(r),v(r,t a )-v'(r,t ) \%) 

= m (r)V[v(r,t Q )-v'(r,t )} . (22) 

Let us assume that Eq. (13) is fulfilled already for k — 0, i.e. that the two 
potentials, v and v' , differ at to- This immediately implies that the derivative 
on the left-hand side of Eq. (22) differs from zero. The two current densities 
j and j' will consequently deviate for t > to. If k is greater than zero the 
equation of motion is applied k + 1 times to yield 



d k+i 



t=t 



n (r)Vufc(r) 



(23) 



The right-hand side of Eq. (23) differs from zero, which again implies that 
j{r,t)^j'{r,t)iovt>t . 

In a second step we prove that j ^ j' implies n ^ n' . To achieve that 
purpose we will make use of the continuity equation 



d_ 

dt 



n(r,t) = -V-j(r,t). 



(24) 



If we write Eq. (24) for the primed and unprimed system and take the dif- 
ference, we arrive at 



d_ 

dt 



(r, t) - n'(r, t)} - -V • [j(r, t) - j'(r, t)] . 



(25) 
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As before, we would like an expression involving the fcth time derivative of 
the external potential. We therefore take the (k + l)st time-derivative of the 
previous equation to obtain (at t = to) 

Qk+2 Qk+1 

OtkT-2 KM) - n'(r,t)] t=t0 = -V • ^ [j(r,t)-j'(r,t)] t=to 

= -V • [no(r)Vu fc (r)] . (26) 

In the last step we made use of Eq. (23). By the hypothesis (13) we have 
Uk(r) ^ const, hence it is clear that if 

V ■ [no(r)V« fc (r)] ^0, (27) 

then n ^ n! , from which follows the Runge-Gross theorem. To show that 
Eq. (27) is indeed fulfilled, we will use the versatile technique of demonstra- 
tion by reductio ad absurdum. Let us assume that V • [no(r)Vufe(r)] = with 
Mfc(r) 7^ const., and look at the integral 

J d 3 r «o(r) [Vw fc (r)] 2 - - J d 3 r u fe (r)V • [n (r)Vu fc (r)] (28) 

n (r)u k (r)Vu k (r) ■ dS . 



This equality was obtained with the help of Green's theorem. The first term 
on the right-hand side is zero by assumption, while the second term vanishes 
if the density and the function Uk(r) decay in a "reasonable" manner when 
r — > oo. This situation is always true for finite systems. We further notice that 
the integrand no(r) [Vttfe(r)] 2 is always positive. These diverse conditions 
can only be satisfied if either the density n or Viife(r) vanish identically 
The first possibility is obviously ruled out, while the second contradicts our 
initial assumption that Uk(r) is not a constant. This concludes the proof of 
the Runge-Gross theorem. 



2.3 Time-dependent Kohn-Sham equations 

As mentioned in section 2.1, the Runge-Gross theorem asserts that all observ- 
ables can be calculated with the knowledge of the one-body density. Nothing 
is however stated on how to calculate that valuable quantity. To circumvent 
the cumbersome task of solving the interacting Schrodingcr equation, Kohn 
and Sham had the idea of utilizing an auxiliary system of non-interacting 
(Kohn-Sham) electrons, subject to an external local potential, «ks[7]. This 
potential is unique, by virtue of the Runge-Gross theorem applied to the non- 
interacting system, and is chosen such that the density of the Kohn-Sham 
electrons is the same as the density of the original interacting system. In the 
time-dependent case, these Kohn-Sham electrons obey the time-dependent 
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Schrodingcr equation 

V 2 



• 9 . , 



, , v KS (r,t) 



<Pi(r,t). (29) 



The density of the interacting system can be obtained from the time-dependent 
Kohn-Sham orbitals 



occ 



n(r,t) = J2 bi( r »*)| 2 • ( 3 °) 

i 

Eq. (29), having the form of a one-particle equation, is fairly easy to solve 
numerically. We stress, however, that the Kohn-Sham equation is not a mean- 
field approximation: If we knew the exact Kohn-Sham potential, t>KSj we 
would obtain from Eq. (29) the exact Kohn-Sham orbitals, and from these 
the correct density of the system. 

The Kohn-Sham potential is conventionally separated in the following 
way 

VK.s(r,t) = V cxt (r,t) + UHartree(r,t) + V xc (r,t) . (31) 

The first term is again the external potential. The Hartree potential accounts 
for the classical electrostatic interaction between the electrons 

vu,rtree(r,t)= [d 3 r'^-4r. (32) 



The last term, the xc potential, comprises all the non-trivial many-body 
effects. In ordinary DFT, v xc is normally written as a functional derivative 
of the xc energy. This follows from a variational derivation of the Kohn- 
Sham equations starting from the total energy. It is not straightforward to 
extend this formulation to the time-dependent case due to a problem related 
to causality [11, 2]. The problem was solved by van Leeuwen in 1998, by using 
the Keldish formalism to defined a new action functional [12], A. The time- 
dependent xc potential can then be written as the functional derivative of 
the xc part of A, 

v KC {r,t) 



5n(r, t) 



(33) 



where r stands for the Keldish pseudo-time. 

Inevitably, the exact expression of w xc as a functional of the density is 
unknown. At this point we are obliged to perform an approximation. It is im- 
portant to stress that this is the only fundamental approximation in TDDFT. 
In contrast to stationary-state DFT, where very good xc functional exist, 
approximations to v xc (r, t) are still in their infancy. The first and simplest of 
these is the adiabatic local density approximation (ALDA), reminiscent of the 
ubiquitous LDA. More recently, several other functional were proposed, from 
which we mention the time-dependent exact-exchange (EXX) functional [13], 
and the attempt by Dobson, Biinner, and Gross[14] to construct an xc func- 
tional with memory. In the following section we will introduce the above 
mentioned functional . 
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2.4 xc functionals 

Adiabatic approximations There is a very simple procedure that allows 
the use of the plethora of existing xc functionals for ground-state DFT in the 
time-dependent theory. Let us assume that v xc [n] is an approximation to the 
ground-state xc density functional. We can write an adiabatic time-dependent 
xc potential as 

w adiabatic (pjt) = 5xc [ n ](r)|„ = „ (t) . (34) 

I.e. we employ the same functional form but evaluated at each time with the 
density n(r,t). The functional thus constructed is obviously local in time. 
This is, of course, a quite dramatic approximation. The functional w xc [n] is 
a ground-state property, so we expect the adiabatic approximation to work 
only in cases where the temporal dependence is small, i.e., when our time- 
dependent system is locally close to equilibrium. Certainly this is not the case 
if we are studying the interaction of strong laser pulses with matter. 

By inserting the LDA functional in Eq. (34) we obtain the so-called adi- 
abatic local density approximation (ALDA) 

C DA (^)-fW|„ = „ M . (35) 

The ALDA assumes that the xc potential at the point r, and time t is equal 
to the xc potential of a (static) homogeneous-electron gas (HEG) of density 
n(r, t). Naturally, the ALDA retains all problems already present in the LDA. 
Of these, we would like to emphasize the erroneous asymptotic behavior of the 
LDA xc potential: For neutral finite systems, the exact xc potential decays 
as — 1/r, whereas the LDA xc potential falls off exponentially Note that 
most of the generalized-gradient approximations (GGAs), or even the newest 
meta-GGAs have asymptotic behaviors similar to the LDA. This problem 
gains particular relevance when calculating ionization yields (the ionization 
potential calculated with the ALDA is always too small), or in situations 
where the electrons are pushed to regions far away from the nuclei (e.g., by 
a strong laser) and feel the incorrect tail of the potential. 

Despite this problem, the ALDA yields remarkably good excitation ener- 
gies (see sections 4.2 and 4.3) and is probably the most used xc functional in 
TDDFT. 



Time-dependent optimized effective potential Unfortunately, when 
one is trying to write v xc as explicit functionals of the density, one encoun- 
ters some difficulties. As an alternative, the so-called orbital-dependent xc 
functionals were introduced several years ago. These functionals are written 
explicitly in terms of the Kohn-Sham orbitals, albeit remaining implicit den- 
sity functionals by virtue of the Runge-Gross theorem. A typical member of 
this family is the exact-exchange (EXX) functional. The EXX action is ob- 
tained by expanding A^ c in powers of e 2 (where e is the electronic charge), 
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and retaining the lowest order term, the exchange term. It is given by the 
Fock integral 



to 



■>»r A = -§2J/ ^/^/^ — i. In ■ (36) 



From such an action functional, one seeks to determine the local Kohn-Sham 
potential through a series of chain rules for functional derivatives. The pro- 
cedure is called the optimized effective potential (OEP) or the optimized 
potential method (OPM) for historical reasons [15, 16]. The derivation of the 
time-dependent version of the OEP equations is very similar to the ground- 
state case. Due to space limitations we will not present the derivation in this 
article. The interested reader is advised to consult the original paper [13], one 
of the more recent publications[17,18], or the chapter by E. Engel contained 
in this volume. The final form of the OEP equation that determines the EXX 
potential is 

OCC ~ti r- 

= E/ dt' d 3 r' MrVj-UxidV)] (37) 
x<p j {r,t)<p*{r',t')G R {rt,r't')+ c.c. 

The kernel, Gr, is defined by 

oo 

iG R (rt, r't') = ]T ^(r, t)ip k {r' , t')6{t - t>) , (38) 
fe=i 

and can be identified with the retarded Green's function of the system. More- 
over, the expression for u x is essentially the functional derivative of the xc 
action in relation to the Kohn-Sham wave-functions 

Uxj (r,t)= J SA *M (39) 

Note that the xc potential is still a local potential, albeit being obtained 
through the solution of an extremely non-local and non-linear integral equa- 
tion. In reality, the solution of Eq. (37) poses a very difficult numerical prob- 
lem. Fortunately, by performing an approximation first proposed by Krieger, 
Li, and Iafratc (KLI) it is possible to simplify the whole procedure, and ob- 
tain an semi-analytic solution of Eq. (37) [19]. The KLI approximation turns 
out to be a very good approximation to the EXX potential. Note that both 
the EXX and the KLI potential have the correct — 1/r asymptotic behavior 
for neutral finite systems. 



A functional with memory There is a very common procedure for the 
construction of approximate xc functional in ordinary DFT. It starts with 
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the derivation of exact properties of w xc , deemed important by physical ar- 
guments. Then an analytical expression for the functional is proposed, such 
that it satisfies those rigorous constraints. We will use this recipe to generate 
a time-dependent xc potential which is non-local in time, i.e. that includes 
the "memory" from previous times [14]. 

A very important condition comes from Galilean invariancc. Let us look 
at a system from the point of view of a moving reference frame whose origin is 
given by x(t). The density seen from this moving frame is simply the density 
of the reference frame, but shifted by x(t) 

n'(r,t) = n(r-x(t),t). (40) 

Galilean invariance then implies [20] 

v xc [n'](r,t) = v xc [n](r - x(t),t) . (41) 

It is obvious that potentials that are both local in space and in time, like the 
ALDA, trivially fulfill this requirement. However, when one tries to deduce 
an xc potential which is non-local in time, one finds condition (41) quite 
difficult to satisfy. 

Another rigorous constraint follows from Ehrcnfcst's theorem which re- 
lates the acceleration to the gradient of the external potential 

^<r> = -<V«ext(r)) • (42) 



For an interacting system, Ehrcnfcst's theorem states 

/ d 3 rrn(r,t) = - J d 3 r n(r,t)Vv cxt (r) . (43) 
In the same way we can write Ehrcnfcst's theorem for the Kohn-Sham system 
Jd 3 rr n(r,t) = - J d 3 r n(r,t)\7v KS (r) . (44) 



dt 2 



By the construction of the Kohn-Sham system, the interacting density is 
equal to the Kohn-Sham density. We can therefore equate the right-hand 
sides of Eq. (43) and (44), and arrive at 

J d 3 r n(r, t)Vv cxt (r) = J d 3 r n{r, t)Vv KS {r, t) . (45) 

If we now insert the definition of the Kohn-Sham potential, Eq. 31, and note 
that / d 3 r n(r, i)VwHartroc('") = 0, we obtain the condition 

jd 3 rn{r 1 t)\7v xc (r 1 t) = Jd 3 rn{r,t)F xc (r,t) = 0, (46) 
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i.e. the total xc force of the system is zero. This condition reflects Newton's 
third law: The xc effects are only due to internal forces, the Coulomb inter- 
action among the electrons, and should not give rise to any net force on the 
system. 

A functional that takes into account these exact constraints can be con- 
structed [14]. The condition (46) is simply ensured by the expression 

F xc (r,t) = -^—jV Jdt'n^{n{r,t'),t-t'). (47) 

The function 7T XC is a pressure-like scalar memory function of two variables. 
In practice, 77 xc is fully determined by requiring it to reproduce the scalar 
linear response of the homogeneous electron gas. Expression (47) is clearly 
non-local in the time-domain but still local in the spatial coordinates. From 
the previous considerations it is clear that it must violate Galilean invariance. 
To correct this problem we use a concept borrowed from hydrodynamics. It 
is assumed that, in the electron liquid, memory resides not with each fixed 
point r, but rather within each separate "fluid element". Thus the element 
which arrives at location r at time t "remembers" what happened to it at 
earlier times t' when it was at locations R(t'\r,t), different from its present 
location r. The trajectory, R, can be determined by demanding that its time 
derivative equals the fluid velocity 

with the boundary condition 

R(t\r,t)=r. (49) 
We then correct the Eq. (47) by evaluating n at point R 

Fxc(r,i) = -r i-r V [dt'n xc {n{R,t'),t-t'). (50) 
n{r, t) J 

Finally, an expression for v xc can be obtained by direct integration of F xc 
(see [14] for details). 



2.5 Numerical considerations 

As mentioned before, the solution of the time-dependent Kohn-Sham equa- 
tions is an initial value problem. At t = to the system is in some initial state 
described by the Kohn-Sham orbitals ifi(r, to)- In most cases the initial state 
will be the ground state of the system (i.e., ipi(r,to) will be the solution of 
the ground-state Kohn-Sham equations). The main task of the computational 
physicist is then to propagate this initial state until some final time, tf. 
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The time-dependent Kohn-Sham equations can be rewritten in the inte- 
gral form 

ifi{r,tf) = U{tf,t )<pi(r,t ) , (51) 



where the time-evolution operator, U, is defined by 



U(t',t) = f exp 



-l J dr H ks (t) 



(52) 



Note that -ffxs i s explicitly time-dependent due to the Hartree and xc po- 
tentials. It is therefore important to retain the time-ordering propagator, T, 
in the definition of the operator U. The exponential in expression (52) is 
clearly too complex to be applied directly, and needs to be approximated 
in some suitable manner. To reduce the error in the propagation from to 
to tf, this large interval is usually split into smaller sub-intervals of length 
A t. The wave-functions are then propagated from to — > to + A t, then from 
to + A t — > to + 2A t and so on. 

The simplest approximation to (56) is a direct expansion of the exponen- 
tial in a power series of A t 



U(t+At,t) w 



\-iH{t + At/2)At 



1=0 



0(At k+L ). 



(53) 



Unfortunately, the expression (53) does not retain one of the most impor- 
tant properties of the Kohn-Sham time-evolution operator: unitarity. In other 
words, if we apply Eq. (53) to a normalized wave-function the result will no 
longer be normalized. This leads to an inherently unstable propagation. Sev- 
eral propagation methods that fulfill the condition of unitarity exist in the 
market. We will briefly mention two of these: a modified Krank-Nicholson 
scheme, and the split-operator method. 



A modified Krank-Nicholson scheme This method is derived by im- 
posing time-reversal symmetry to an approximate time-evolution operator. 
It is clear that we can obtain the state at time t + At/2 either by forward 
propagating the state at t by A t/2, or by backward propagating the state at 
t + At 

<p(t + A t/2) = U(t + A t/2, t)ip(t) 

= U(t- At/2,t+ At)tp(t + At) . (54) 

This equality leads to 

ip{t + At) = U(t + A t/2, t + A t)U(t + A t/2, t)<p(t) , (55) 

where we used the fact that the inverse of the time-evolution operator U~ 1 (t+ 
At, t) = U (t — At, t). To propagate a state from t to t+At we follow the steps: 
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i) Obtain an estimate of the Kohn-Sham wave-functions at time t + A t by 
propagating from time t using a "low quality" formula for U(t + At,t). The 
expression (53) expanded to third or forth order is well suited for this purpose. 

ii) With these wave-functions construct an approximation to H(t + At) and 
to U(t + A t/2, t + At), iii) Apply Eq. (55). This procedure leads to a very 
stable propagation. 



The split-operator method In a first step we neglect the time-ordering in 
Eq. (52), and approximate the integral in the exponent by a trapezoidal rule 



U(t + A t, t) w exp -iH KS (t)At 



exp 



-i(T+V KS )At 



(56) 



We note that the operators exp (— iVksA t) and exp (— if A t) are diagonal 

respectively in real and Fourier spaces, and therefore trivial to apply in those 
spaces. It is possible to decompose the exponential (56) into a form involving 
only these two operators. The two lowest order decompositions are 



exp 
and 
exp 



-iCf+Vks)^* = exp (-if At) exp (-iV KS At) +0(At 2 ) (57) 

-i(T + V K s)A i] = exp (-jf^J exp (-iVks^ t) exp (-i^4r 

+0(At 3 ). (58) 



For example, to apply the splitting (58) to ip(r, t) we start by Fourier trans- 
forming the wave-function to Fourier space. We then apply exp (— ifA±) to 
ip(k,t) and Fourier transform back the result to real space. We proceed by 
applying exp (-iVA t) , Fourier transforming, etc. This method can be made 
very efficient by the use of fast-Fourier transforms. 

As a better approximation to the propagator (52) we can use a mid-point 
rule to estimate the integral in the exponential 



U(t + A t, t) w exp -iH K s(t + A t/2) At 



(59) 



It can be shown that the same procedure described above can be applied 
with only a slight modification: The Kohn-Sham potential has to be updated 
after applying the first kinetic operator [21]. 



3 Linear response theory 
3.1 Basic theory 

In circumstances where the external time-dependent potential is small, it 
may not be necessary to solve the full time-dependent Kohn-Sham equations. 
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Instead perturbation theory may prove sufficient to determine the behavior 
of the system. We will focus on the linear change of the density, that allows 
us to calculate, e.g., the optical absorption spectrum. 

Let us assume that for t < t the time-dependent potential vtd is zero - 
i.e. the system is subject only to the nuclear potential, - and furthermore 
that the system is in its ground-state with ground-state density n^K At to 
we turn on the perturbation, so that the total external potential now 
consists of w ox t = + . Clearly will induce a change in the density. 
If the perturbing potential is sufficiently well-behaved (like almost always in 
physics), we can expand the density in a perturbative series 

n(r, t) = n (0) (r) + n (1) (r,t) + n {2) (r, t) + ■ ■ ■ (60) 

where is the component of n(r,t) that depends linearly on 
depends quadratically, etc. As the perturbation is weak, we will only be con- 
cerned with the linear term, In frequency space it reads 



nW(r,^) = J "d 3 r' ' X (r,r' \uj)v^(r' » . (61) 

The quantity \ is the linear density-density response function of the sys- 
tem. In other branches of physics it has other names, e.g., in the context of 
many-body perturbation theory it is called the reducible polarization func- 
tion. Unfortunately, the evaluation of x through perturbation theory is a 
very demanding task. We can, however, make use of TDDFT to simplify this 
process. 

We recall that in the time-dependent Kohn-Sham framework, the density 
of the interacting system of electrons is obtained from a fictitious system of 
non-interacting electrons. Clearly, we can also calculate the linear change of 
density using the Kohn-Sham system 



nW(r,w) 



ydVxKs(r,r>)«W(r>). (62) 



Note that the response function that enters Eq. (62), xks> is the density re- 
sponse function of a system of non-interacting electrons and is, consequently, 
much easier to calculate than the full interacting \- I n terms of the unper- 
turbed stationary Kohn-Sham orbitals it reads 

* K3(r ' r ' W) = ^+ ^ {h ^ oj-te-eO + ii? ' (63) 

where f m is the occupation number of the m orbital in the Kohn-Sham 
ground-state. Note that the Kohn-Sham potential, t>KS, includes all powers 
of the external perturbation due to its non-linear dependence on the density. 
The potential that enters Eq. (62) is however just the linear change of vks, 
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v£s- This latter quantity can be calculated explicitly from the definition of 
the Kohn-Sham potential 

«^(r,t) = v^(r,t) + v^ Itiee (r,t) + v£(r,t) ■ (64) 

The variation of the external potential is simply , while the change in the 
Hartree potential is 



[artrec( r > ^ 



•J 1 

J Hart 



r — r'\ 



(65) 



Finally v x ^(r,t) is the linear part in of the functional u xc [n], 



«£ ) (r,*) = f#f*'j$$» ll) ir'J). (66) 



It is useful to introduce the exchange-correlation kernel, / xc , defined by 

<5v xc (r,t) 



/xc(rt,rY) 



Sn(r',t') 



(67) 



The kernel is a well know quantity that appears in several branches of the- 
oretical physics. E.g., evaluated for the electron gas, / xc is, up to a factor, 
the "local-field correction" ; To emphasize the correspondence to the effective 
interaction of Landau's Fermi-liquid theory, / xc plus the bare Coulomb in- 
teraction is sometimes called the "effective interaction" , while in the theory 
of classical liquids the same quantity is referred to as the Ornstein-Zernicke 
function. 

Combining the previous results, and transforming to frequency space we 
arrive at: 



(r,u) = j d 3 r' XKs(r,r' ,u)v {1 \r' ,w) 



+ / d 6 x / dVxKs(r,x,w) 



7l + f xc {x,r',uj) 



\x — r 

From Eq. (61) and Eq. (68) trivially follows the relation 
X(r,r',u) = XKs{r,r',u) + 



d x Id x x(r, x, uj) 



IX — X 



7j + f*c{x,x',u)) 



(68) 
n«(r>). 

(69) 

XKs{x',r',u) . 



This equation is a formally exact representation of the linear density response 
in the sense that, if we possessed the exact Kohn-Sham potential (so that we 
could extract / xc ), a self-consistent solution of (69) would yield the response 
function, x, of the interacting system. 
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3.2 The xc kernel 

As we have seen in the previous section, the main ingredient in linear response 
theory is the xc kernel. / xc , as expected, is a very complex quantity that 
includes - or, in other words, hides - all non-trivial many-body effects. Many 
approximate xc kernels have been proposed in the literature over the past 
years. The most ancient, and certainly the simplest is the ALDA kernel 

f^{rty t>) = 8{r r')S(t t>) / x H c EG ML„ (T ,, t) , (70) 

where 

/x H c EG W ^ ^x H c EG (") (71) 

is just the derivative of the xc potential of the homogeneous electron gas. 
The ALDA kernel is local both in the space and time coordinates. 

Another commonly used xc kernel was derived by Pctcrsilka et al. in 1996, 
and is nowadays referred to as the PGG kernel [22]. Its derivation starts from 
a simple analytic approximation to the EXX potential. This approximation, 
that is called the Slater approximation in the context of Hartree-Fock theory, 
only retains the leading term in the expression for EXX, which reads 

occ i , n,2 

^ GG MHE^7T [Ux*(r, t)+cx] . (72) 

Using the definition (67) and after some algebra we arrive at the final form 
of the PGG kernel 

f PGG (rt r't') - -Sit - 1') 1 1 KrVfcfrMfrOI 2 (J3) 
; x (rt,rt)- d(t *; 2 | p _ r ,| n ( P ) n ( P /) ' ^ 

As in the case of the ALDA, the PGG kernel is local in time . 

Noticing the crudeness of the ALDA, especially the complete neglect of 
any frequency dependence, one could expect it to yield very inaccurate re- 
sults in most situations. Surprisingly, this is not the case as we will show in 
section 4. To understand this numerical evidence, we have to take a step back 
and study more thoroughly the properties of the xc kernel for the homoge- 
neous electron gas[l]. 

In this simple system /™ G only depends on r — r' and on t — t' , so 
it is convenient to work in Fourier space. Our knowledge of the function 
f^ G (q,oj) is quite limited. Several of its exact features can nevertheless 
be obtained through analytical manipulations. The zero frequency and zero 
momentum limit is given by 

lim /x HEG (g, co = 0) = ^ K c EG (n)] ee /„(„) , (74) 
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where e xc EG , the xc energy per particle of the homogeneous electron gas, is 
known exactly from Monte-Carlo calculations [23]. Also the infinite frequency 
limit can be written as a simple expression 



lim f™- G (q,uj = oo) = -^n§ 4~ 
= fooin) . 



HEGf \1 d r HEG/ y 



,2/3 



From these two expression, one can prove that the zero frequency limit is 
always smaller than the infinite frequency limit, and that both these quanti- 
ties are smaller than zero (according to the best approximations known for 
£ x H c EG ), i.e. 

fo(n) < fooin) < (76) 

From the fact that f xc is a real function when written in real space and in 
real time one can deduce the following symmetry relations 



$tf™ G (q,co) = $lf™ G ( q ,-u,) (77) 



From causality follow the Kramers-Kronig relations: 
Rf x H c EG (<7, ") - / x H c EG ( 9 , oo) = P T ^ 5/ " EG(g '" ) (78) 

J-00 I 1 UJ - OJ' 

where V denotes the principal value of the integral. Note that, as the infinite 
frequency limit of the xc kernel is different from zero, one has to subtract 
,/™ G (q, 00) in order to apply the Kramers-Kronig relations. 

Furthermore, by performing a perturbative expansion of the irreducible 
polarization to second order in e 2 , one finds 

Jim 9/^ = 0,0,) = -^. (79) 

The real part can be obtained with the help of the Kramers-Kronig relations 

Jim Rf x H c EG (<7 = 0, cj) = /«,(„) + ^-3^ . (80) 

It is possible to write a analytical form for the long-wavelength limit of the 
imaginary part of / xc that incorporates all these exact limits [24] 

S /x HEG ( g = 0,^ a{n)u „ . (81) 
(l + /3(n)w 2 )f 




Fig. 1. Real and imaginary part of the parametrization for /™ . Figure reproduced 
from Rcf. [25]. 



The coefficients a and [3 are functions of the density, and can be determined 
uniquely by the zero and high frequency limits. A simple calculation yields 

a(n) = -4[/oo(n)-/o(n)]* (82) 
P(n) = B[f oo (n)-f (n)]i, (83) 

where A, B > and independent of n. Once again, by applying the Kramers- 
Kronig relations we can obtain the corresponding real part of /™ G 



(^) 



1 + r n f l-r J_\ _ 1 -r n ( \ r ! 



V2J 2 V 2 'V2 



where r = + /3to 2 and S and II are the elliptic integrals of second and 
third kind. In Fig. 1 we plot the real and imaginary part of /™ G for two 
different densities (r s = 2 and r s = 4, where r s is the Wigner-Seitz radius, 
1/n = 47rrg/3). The ALDA corresponds to approximating these curves by 
their zero frequency value. For very low frequencies, the ALDA is naturally a 
good approximation, but at higher frequencies it completely fails to reproduce 
the behavior of /™ G . 

To understand how the ALDA can yield such good excitation energies, 
albeit exhibiting such a mediocre frequency dependence, we will look at a spe- 
cific example, the process of photo-absorption by an atom. At low excitation 
frequencies, we expect the ALDA to work. As we increase the laser frequency, 
we start exciting deeper levels, promoting electrons from the inner shells of 
the atom to unoccupied states. The atomic density increases monotonically 
as we approach the nucleus. / xc corresponding to that larger density (lower 
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r s ) has a much weaker frequency dependence, and is much better approxi- 
mated by the ALDA than the low density curve. In short, by noticing that 
high frequencies are normally related to high densities, we realize that for 
practical applications the ALDA is often a reasonably good approximation. 
One should however keep in mind that these are simple heuristic arguments 
that may not hold in a real physical system. 

4 Excitation energies 

4.1 DFT techniques to calculate excitations 

In this section we will present a short overview of the several techniques to 
calculate excitation energies that have appeared in the context of DFT over 
the past years. Indeed quite a lot of different approaches have been tried. 
Some are more or less ad hoc, others rely on a solid theoretical basis. More- 
over, the degree of success varies considerably among the different techniques. 
The most successful of all is certainly TDDFT that has become the de facto 
standard for the calculation of excitations for finite systems. We will leave 
the discussion of excitation energies in TDDFT to the following sections, and 
concentrate for now on the "competitors" . The first group of methods is based 
on a single determinant calculation, i.e. only one ground-state like calcula- 
tion is performed, subject to the restriction that the Kohn-Sham occupation 
numbers are either or 1. 

As a first approximation to the excitation energies, one can simply take 
the differences between the ground-state Kohn-Sham eigenvalues. This pro- 
cedure, although not entirely justifiable, is often used to get a rough idea of 
the excitation spectrum. We stress that the Kohn-Sham eigenvalues (as well 
as the Kohn-Sham wave-functions) do not have any physical interpretation. 
The exception is the eigenvalue of the highest occupied state that is equal to 
minus the ionization potential of the system [26]. 

The second scheme is based on the observation that the Hohenberg-Kohn 
theorem and the Kohn-Sham scheme can be formulated for the lowest state 
of each symmetry class [27]. In fact, the single modification to the standard 
proofs is to restrict the variational principle to wave-functions of a specific 
symmetry. The unrestricted variation will clearly yield the ground-state. The 
states belonging to different symmetry classes will correspond to excited 
states. The excitations can then be calculated by simple total-energy differ- 
ences. This approach suffers from two serious drawbacks: i) Only the lowest 
lying excitation for each symmetry class is obtainable, ii) The xc functional 
that now enters the Kohn-Sham equations depends on the particular sym- 
metry we have chosen. As specific approximations for a symmetry dependent 
xc functional are not available, one is relegated to use ground-state function- 
al. Unfortunately the excitation energies calculated in this way are only of 
moderate quality. 
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Another promising method was recently proposed by A. G6rling[28]. The 
so-called generalized adiabatic connection Kohn-Sham formalism is no longer 
based on the Hohcnbcrg-Kohn theorem but on generalized adiabatic connec- 
tions associating a Kohn-Sham state with each state of the real system. This 
formalism was later extended to allow for the proper treatment of symme- 
try of the Kohn-Sham states[29]. The quality of the results obtained so far 
with this procedure varies: For alkali atoms the agreement with experimental 
excitation energies is quite good [28], but for the Carbon atom and the CO 
molecule the situation is considerably worse [29]. We note however that this 
method is still in its infancy so further developments can be expected in the 
near future. 

It is also possible to calculate excitation energies from the ground-state 
energy functional. In fact, it was proved by Perdew and Levy [30] that "every 
extremum density rii (r) of the ground-state energy functional E v [n] yields the 
energy Ei of a stationary state of the system." The problem is that not every 
excited-state density, rij(r), corresponds to an extremum of i£„[n], which 
implies that not all excitation energies can be obtained from this procedure. 

The last member of the first group of methods was proposed by Ziegler, 
Rauk and Baerends in 1977[31] and is based on an idea borrowed from multi- 
configuration Hartree-Fock. The procedure starts with the construction of 
many-particle states with good symmetry, tf'i, by taking a finite superposition 
of states 

a 

where <P a are Slater determinants of Kohn-Sham orbitals, and the coefficients 
Ci a are determined from group theory. Through a simple matrix inversion we 
can express the determinants as linear combinations of the many-body wave- 
functions 

*/J = X>«^ - - (86) 

3 

By taking the expectation value of the Hamiltonian in the state $p we arrive 
at 

($ \H\<P fj )=J2M 2 E 3 , (87) 

3 

where Ej is the energy of the many-body state The "recipe" to calculate 
excitation energies is then: a) Build <Pp from n Kohn-Sham orbitals (not 
necessarily the lowest); b) Make an ordinary Kohn-Sham calculation for each 
$p, and associate the corresponding total energy E^ FT with H 
c) Determine Ej by solving the system of linear equations, Eq. (87). 

This methods works quite well in practice, and was frequently used in 
quantum chemistry till the advent of TDDFT. We should nevertheless indi- 
cate two of its limitations: i) The decomposition (85) is not unique and the 
system of linear equations can be under- or over determined, ii) The whole 
procedure of the "recipe" is not rigorously founded. 
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The next technique, the so-called ensemble DFT, makes use of frac- 
tional occupation numbers. Ensemble DFT, first proposed by Thcophilou 
in 1979[32], evolves around the concept of an ensemble. In the simplest case 
it consists of a "mixture" of the ground state, &i, and the first excited state, 
<?2, described by the density matrix[33-35] 

Z?=(l-w)|#l)<!Pl|+W |^ 2 ) ON , (88) 

where the weight, uj, is between and 1/2 (in this last case the ensemble 
is called "equiensemble" ) . We can further define the ensemble energy and 
density 

E(u>) = (1 - uS)E 1 + w E 2 (89) 
n w (r) = (l-w)m(r)+wn 2 (r). (90) 

At oj — the ensemble energy clearly reduces to the ground-state energy. 
Using the ensemble density it is possible to construct a DFT, i.e. to prove a 
Hohenbcrg-Kohn theorem and construct a Kohn-Sham scheme. The main fea- 
tures of the Kohn-Sham scheme arc: i) The one-body orbitals have fractional 
occupations determined by the weight u>. ii) The xc functional depends on 
the weight, E xc (u). To calculate the excitation energies from ensemble DFT 
we can follow two paths. The first involves obtaining the ground-state energy 
and the ensemble energy for some fixed u>, from which the excitation energy 
Ei — E 2 trivially follows 

E 2 -Ei = E ^~ E ^. (91) 

The second path is obtained by taking the derivative of Eq. (89) 

dE(uj) 
duj 

It is then possible to prove 



= E 2 - Ei . (92) 



' • ■ ' • 1 — e„, — e,„ - - 



duo 



(93) 



Naturally, we need approximations to the xc energy functional, E xc (oj). An 
ensemble LDA was developed for the equiensemble by W. Kohn in 1986 [36], 
by treating the ensemble as a reminiscent of a thermal ensemble. He then re- 
lated E xc (oj) to the finite temperature xc energy of the homogeneous electron 
gas by equating the entropies of both systems. Unfortunately, the results ob- 
tained with this functional were not very encouraging. A promising approach, 
recently proposed, is the use of orbital functionals within an ensemble OEP 
method[37,38]. 
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4.2 Full solution of the Kohn-Sham equations 

One of the most important uses of TDDFT is the calculation of photo- 
absorption spectra. This problem can be solved in TDDFT either by prop- 
agating the time-dependent Kohn-Sham equations [39] or by using linear- 
response theory. In this section we will be concerned by the former, relegating 
the latter to the next section. 

Let <pj (r) be the ground-state Kohn-Sham wave-functions for the system 
under study. We prepare the initial state for the time-dependent propagation 
by exciting the electrons with the electric field v(r,t) = —kox„S(t), where 
x v = x,y,z. The amplitude ko must be small in order to keep the response 
of the system linear and dipolar. Through this prescription all frequencies of 
the system are excited with equal weight. At t = + the initial state for the 
time evolution reads 



The Kohn-Sham orbitals are then further propagated during a finite time. 
The dynamical polarizability can be obtained from 



In the last expression 5n(r,u)) stands for the Fourier transform of n(r,t) — 
n(r), where h{r) is the ground-state density of the system. The quantity that 
is usually measured in experiments, the photo-absorption cross-section, is 
essentially proportional to the imaginary part of the dynamical polarizability 
averaged over the three spatial directions 



where c stands for the velocity of light. Although computationally more de- 
manding than linear-response theory, this method is very flexible, and is 
easily extended to incorporate temperature effects, non-linear phenomena, 
etc. Note also that this approach only requires an approximation to the xc 
potential. 

To illustrate the method, we present, in Fig. 2, the excitation spectrum 
of benzene calculated within the LDA/ALDA 3 . The agreement with experi- 
ment is quite remarkable, especially when looking at the 7r — > 7r* resonance at 

3 We will use the notation "A/B" consistently throughout the rest of this article to 
indicate that the ground-state xc potential used to calculate the initial state was 
"A", and that this state was propagated with the time-dependent xc potential 
"B" . In the case of linear-response theory, "B" will denote the xc kernel. 




cxp [ik x v ] <f>j(r) . 



(94) 




(95) 




(96) 
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Fig. 2. Optical absorption of the benzene molecule. Experimental results from 
Ref. [40]. Figure reproduced from Ref. [41]. 

around 7eV. The spurious peaks that appear in the calculation at higher en- 
ergies are artifacts caused by an insufficient treatment of the unbound states. 
We furthermore observe that such good results are routinely obtained when 
applying the LDA/ALDA to several finite systems, from small molecules to 
metallic clusters and biological systems. 

4.3 Excitations from linear-response theory 

The first self-consistent solution of the linear response equation (69) was 
performed by Zangwill and Soven in 1980 using the LDA/ALDA [42]. Their 
results for the photo-absorption spectrum of Helium for energies just above 
the ionization threshold are shown in Fig. 3. Once more the theoretical curve 
compares very well to experiments. 

Unfortunately, a full solution of Eq. (69) is still quite difficult numerically. 
Besides the large effort required to solve the integral equation, we need as 
an input the non-interacting response function. To obtain this quantity it is 
usually necessary to perform a summation over all states, both occupied and 
unoccupied [cf. Eq. (63)]. Such summations are sometimes slowly convergent 
and require the inclusion of many unoccupied states. There are however ap- 
proximate frameworks that circumvent the solution of Eq. (69) . The one we 
will present in the following was proposed by Petersilka et al. [22]. 
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Fig. 3. Total photo- absorption cross-section of Xenon versus photon energy in the 
vicinity of the 4d threshold. The solid line represents TDDFT calculations and the 
crosses are the experimental results of Ref. [43]. Figure adapted from Ref. [42]. 



The density response function can be written in the Lehmann represen- 
tation 



X(r,r',u) = lim V 
?j^o+ — 



(0\p(r) \m) (m\p(r') |0 
u - (E m - E ) + ir} 



\p(r')\m) (m\p(r)\0) 



oj + (E m - Eq) + if] 



(97) 

where \m) is a complete set of many-body states with energies E m . From this 
expansion it is clear that the full response function has poles at frequencies 
that correspond to the excitation energies of the interacting system 



ft — E m — En . 



(98) 



As the external potential does not have any special pole structure as a func- 
tion of u>, Eq. (61) implies that also n^(r,uj) has poles at the excitation 
energies, Q. On the other hand, %ks has poles at the excitation energies of 
the non-interacting system, i.e. at the Kohn-Sham orbital energy differences 
ej -e k [cf. Eq. (63)]. 

By rearranging the terms in Eq. (68) we obtain the fairly suggestive equa- 
tion 

J d 3 r' [S(r - r') - S(r,r',w)] n (1) (r» = J d 3 r' XKS (r, r', Lo)v {1 \r\ u>) 

(99) 

where the function r; is defined by 



S (r, r',u>) = J d 3 r" xks^, r",u>) 



hc(r",r',Lj) 



(100) 
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As noted previously, in the limit lu — > Q the linear density rS 1 ^ has a pole, 
while the right-hand side of Eq. (100) remains finite. For the equality (100) 
to hold, it is therefore required that E has zero eigenvalues at the excitation 
energies Q, i.e. A(w) — > 1 when u) — ► f2, where X(ui) is the solution of the 
eigenvalue equation 



d 3 r' E(r, r',w)£(r',w) = \(w)£(r,w) . 



(101) 



This is a rigorous statement, that allows the determination of the excitation 
energies of the systems from the knowledge of xks an d / X c- It is possible 
to transform this equation into another eigenvalue equation having the true 
excitation energies of the system, fi, as eigenvalues [44]. We start by defining 
the quantity 



OfcM = J d 3 r'J d 3 r" y*{r")y k {r") 



1 



- + U c (r",r',cu) 



With the help of Qk Eq. (101) can be rewritten in the form 



E 

jk 



{fa - fj)<Pj(r)(p* k (r) 

U - ( e 3 - e k) + ill 



€(r',w). 
(102) 

(103) 



By solving this equation for £(r,w) and inserting the result into Eq. (102), 
we arrive at 



(104) 



where we have defined the matrix clement 



M jk j tk ,(u) = (fk> - fr) J d 3 r J dV^(r)^(r)^(r'K(r') 



\r — r i 

Introducing the new eigenvector 

0jk = 



7i + fa c {r,r',w) 



CM 



Q - (eji - e k >) 



(105) 



(106) 



taking the r\ — > limit, and by using the condition A(J?) = 1, it is straight- 
forward to recast Eq. (104) into the eigenvalue equation 

^ [Sjj'Skk'i^j' - efc') + Mjfc ) j/fc/(/2)] /3j/fe/ = flfijk ■ (107) 
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It is also possible to derive an operator whose eigenvalues are the square 
of the true excitation energies, thereby reducing the dimension of the ma- 
trix equation (107) [45]. The oscillator strengths can be obtained from the 
eigenfunctions of the operator. 

The eigenvalue equation, Eq. (107), can be solved in several different ways. 
For example, it is possible to expand all quantities in a suitable basis and 
solve numerically the resulting matrix-eigenvalue equation. As an alternative, 
we can perform a Laurent expansion of the response function around the 
excitation energy 

( i \ v foo{r)v* {r')ip ko {r')ip* k {r) 
XKs( r , r = lim h higher order . (108) 

»7-»0+ oj - (t jo - t ka ) + irj 

By neglecting the higher-order terms, a simple manipulation of Eq. (101) 
yields the so-called single-pole approximation (SPA) to the excitation energies 

Q = At + K(Ae) , (109) 

where At is the difference between the Kohn-Sham eigenvalue of the unoc- 
cupied orbital jo and the occupied orbital fco, 

At = t jo -t ko (110) 

and if is a correction given by 

K(At) = 2^J d 3 r J d 3 r' ^ (r)ipl(r')<p ko (r'M (r) x (111) 

l^-L_ + / xc(r y,z}e) • 

Although not as precise as the direct solution of the eigenvalue equation, 
Eq. (107), this formula provides us with a simple and fast way to calculate 
the excitation energies. 

To assert how well this approach works in practice we list, in Tab. 1, the 
: S — > P excitation energies for several atoms [22]. Surprisingly perhaps, the 
eigenvalue differences, At, are already of the proper order of magnitude. For 
other systems they can be even much closer (cf. Tab. 3). Adding the correc- 
tion K then brings the numbers indeed very close to experiments for both xc 
functionals tried. We furthermore notice that the EXX/PGG functional gives 
clearly superior results than the LDA/ALDA. This is related to the differ- 
ent quality of the unoccupied states generated with the two ground-state xc 
functionals. The unoccupied states typically probe the farthest regions from 
the system, where the LDA potential exhibits severe deficiencies (as previ- 
ously mentioned in section 2.4). As the EXX potential does not suffer from 
this problem, it yields better unoccupied orbitals and consequently better 
excitation energies. 
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Atom 


ZifLDA ^LDA/ALDA 


Z\eexx 


^EXX/PGG 


i?exp 


Be 


0.129 


0.200 


0.130 


0.196 


0.194 


Mg 


0.125 


0.176 


0.117 


0.164 


0.160 


Ca 


0.088 


0.132 


0.079 


0.117 


0.108 


Zn 


0.176 


0.239 


0.157 


0.211 


0.213 


Sr 


0.082 


0.121 


0.071 


0.105 


0.099 


Cd 


0.152 


0.214 


0.135 


0.188 


0.199 



Table 1. X S — > 1 P excitation energies for selected atoms. J? exp denotes the experi- 
mental results from [46]. All energies are in Hartrees. Table adapted from Ref. [22]. 

In Tab. 2 we show the excitation energies of the CO molecule. This case is 
slightly more complicated than the previous example due to the existence of 
degeneracies in the eigenspectrum of the CO molecule. Although the Kohn- 
Sham eigenvalue differences are equal for all transitions involving degenerate 
states, the true excitation energies depend on the symmetry of the initial and 
final many-body states. As is clearly seen from the table, this splitting of the 
excitations is correctly described by the correction factor, K. 



State 






(? SPA 

LDA/ALDA 


pfull 

"LDA/ALDA 


exp. 


A L n 5a 




2tt 0.2523 


0.3268 


0.3102 


0.3127 


a 3 n 






0.2238 


0.2214 


0.2323 


B L S+ 5a 




6a 0.3332 


0.3389 


0.3380 


0.3962 


b 3 E+ 






0.3315 


0.3316 


0.3822 


i L z- 






0.3626 


0.3626 


0.3631 








0.3626 


0.3626 


0.3631 


a' 3 E+ Itt 




2tt 0.3626 


0.3181 


0.3149 


0.3127 


D M 






0.3812 


0.3807 


0.3759 


d 3 4 






0.3404 


0.3396 


0.3440 


c 3 n 4a 




2tt 0.4388 


0.4204 


0.4202 


0.4245 


E l n Itt 




6a 0.4436 


0.4435 


0.4435 


0.4237 



Table 2. Excitation energies for the CO molecule, ^lda/alda are the LDA/ALDA 
excitation energies obtained from Eq. (109), and ^lda/alda are obtained from 
the solution of Eq. (107) neglecting continuum states, "exp" are the experimental 
results from Ref. [47]. All energies are in Hartrees. Table reproduced from Ref. [48]. 



We remember that several approximations have been made to produce the 
previous results. First a static Kohn-Sham calculation was performed with 
an approximate v KC . Then the resulting cigcnfunctions and eigenvalues were 
used in Eq. (109) to obtain the excitation energies. In the last step, we used 
an approximate form for the xc kernel, f xc , and we neglected the higher order 
terms in the Laurent expansion of the response functions. To assert which of 
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these approximations is more important, we can look at the lowest excitation 
energies of the He atom. For this simple system the exact stationary Kohn- 
Sham potential is known [49], so we can eliminate the first source of error. We 
can then test different approximations for / xc , both by performing the single- 
pole approximation or not. The results are summarized in Tab. 3. We first 
note that the quality of the results is almost insensitive to the xc kernel used. 
Both using the ALDA or the PGG yield the same mean error. This statement 
seems to hold not only for atoms but also to molecular systems [50]. From the 
table it is also clear that the SPA is an excellent approximation and that the 
calculated excitation energies are in very close agreement to the exact values. 
Why, and under which circumstances this is the case is discussed in detail 
in Refs. [51,52]. This leads us to conclude that the crucial approximation to 
obtain excitation energies in TDDFT is the choice of the static xc potential 
used to calculate the Kohn-Sham eigcnfunctions and eigenvalues. 

exact/ALDA (xc) exact/PGG 



State ko 




jo Ae KS SPA 


full 


SPA full exact 


2 3 S 


Is 




2s 0.7460 0.7357 


0.7351 


0.7232 0.7207 0.7285 


2 X S 






0.7718 


0.7678 


0.7687 0.7659 0.7578 


3 3 S 


Is 




3s 0.8392 0.8366 


0.8368 


0.8337 0.8343 0.8350 


3 X S 






0.8458 


0.8461 


0.8448 0.8450 0.8425 


4 a S 


Is 




4s 0.8688 0.8678 


0.8679 


0.8667 0.8671 0.8672 


4 X S 






0.8714 


0.8719 


0.8710 0.8713 0.8701 


2 3 P 


Is 




2p 0.7772 0.7702 


0.7698 


0.7693 0.7688 0.7706 


2 J P 






0.7764 


0.7764 


0.7850 0.7844 0.7799 


3 3 P 


Is 




3s 0.8476 0.8456 


0.8457 


0.8453 0.8453 0.8456 


3 X P 






0.8483 


0.8483 


0.8500 0.8501 0.8486 


4 3 P 


Is 




4s 0.8722 0.8714 


0.8715 


0.8712 0.8713 0.8714 


4 J P 






0.8726 


0.8726 


0.8732 0.8733 0.8727 


Mean abs. dev. 0.0011 
Mean % error 0.15% 


0.0010 
0.13% 


0.0010 0.0010 
0.13% 0.13% 



Table 3. Comparison of the excitation energies of neutral Helium, calculated from 
the exact xc potential[49] by using approximate xc kernels. SPA stands for "sin- 
gle pole approximations", while "full" means the complete solution of Eq. (107). 
The exact values are from a non-relativistic variational calculation [53]. The mean 
absolute deviation and mean percentage errors also include the transitions from 
the Is until the 9s and 9p states. All energies are in Hartrees. Table adapted from 
Ref. [17]. 



4.4 When does it not work? 

In the previous sections we showed the results of several TDDFT calcula- 
tions, most of them agreeing quite well with experiment. Clearly no physical 



Time-dependent Density Functional Theory 31 

theory works for all systems and situations, and TDDFT is not an exception. 
It is the purpose of this section to show some examples where the theory does 
not work. However, before proceeding with our task, we should specify what 
we mean by "failures of TDDFT" . TDDFT is an exact reformulation of the 
time-dependent many-body Schrodingcr equation - it can only fail in situa- 
tions where quantum-mechanics also fails. The key approximation made in 
practical applications is the approximation for the xc potential. Errors in the 
calculations should therefore be imputed to the functional used. As a large 
majority of TDDFT calculations use the ALDA or the adiabatic GGA, we will 
be mainly interested in the errors caused by these approximate functional. 
Furthermore, and as we already mentioned in the previous section, there are 
usually two sources for the errors in the calculation: i) The functional used 
to obtain the Kohn-Sham ground-state; ii) The approximate time-dependent 
xc potential. In any discussion on the errors of TDDFT the effects of these 
two sources have to be clearly separated. With these arguments in mind let 
us then proceed. 

Our first example is the calculation of optical properties of long conju- 
gated molecular chains [54]. Using the local or gradient-corrected approxima- 
tions can give overestimations of several orders of magnitude. The problem 
is related to a non-local dependence of the xc potential: In a system with 
an applied electric field, the exact xc potential develops a linear part that 
counteracts the applied field[54,55]. This term is completely absent in both 
the LDA and the GGA, but is present in more non-local functionals like the 
EXX. 

A related problem occurs in solids [56]. In fact, TDLDA does not work 
properly for the calculation of excitations of non-metallic solids, especially in 
systems like wide-band gap semiconductors. For infinite systems, the Coulomb 
potential is (in momentum space) Air/q 2 . It is then clear from the response 
equation (69) that if / xc is to correct the non-interacting response for g 1 — > it 
will have to contain a term that behaves asymptotically as l/q 2 when q — > 0. 
This is not the case for the local or gradient-corrected approximations. Sev- 
eral attempts have been done to correct this problem from which we mention 
Refs. [57-60]. 

As other faults of the TDLDA we can mention the incorrect excitation 
energies of the stretched H2 molcculc[61], the large error in the calculation 
of singlet-triplet separation energies [62], the underestimate of the onset of 
absorption for some clusters [50], etc. Despite these failures, we would like to 
emphasize that TDLDA does work very well for the calculation of excitations 
in a large class of systems. 
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5 Atoms and molecules in strong laser fields 
5.1 What is a "strong" laser? 

Before discussing the behavior of atoms and molecules in strong laser fields, 
we have to specify what the adjective "strong" means in this context. The 
electric field that an electron feels in a Hydrogen atom, at the distance of one 
Bohr from the nucleus, is 

E= — !— =5.1 x 10 9 V/m. (112) 
47re ao 

The laser intensity that corresponds to this field is given by 

/=ie c£ 2 = 3.51xl0 16 W/cm 2 . (113) 

We can clearly consider a laser to be "strong" when its intensity becomes 
comparable to (113). In this regime, perturbation theory is no longer ap- 
plicable, and the theorist has to resort to non-perturbative methods. When 
approaching these high intensities, a wealth of non-linear phenomena appear, 
like multi-photon ionization, above threshold ionization (ATI) , high harmonic 
generation, etc. 

The fact that allowed systematic investigation of these high-intensity phe- 
nomena was the remarkable evolution in laser technology during the past four 
decades. Through a series of technological breakthroughs, scientists were able 
to boost the peak intensity of pulsed lasers from 10 9 W/cm 2 in the 1960s, to 
more than 10 21 W/cm 2 of the current systems - 12 orders of magnitude! Be- 
sides this increase in laser intensity, very short pulses - sometimes of the order 
of hundreds of attoseconds (las = 10~ 18 s) - became available at ultravio- 
let or soft X-ray frequencies [63,64]. In the present context we are concerned 
mainly with intensities in the range 10 13 — 10 16 W/cm 2 . For higher intensities 
many-body effects associated with the electron-electron interaction - which 
are the main interest of DFT - become less and less important due to the 
strongly dominant external field. 

TDDFT is a tool particularly suited for the study of systems under the 
influence of strong lasers. We recall that the time-dependent Kohn-Sham 
equations yield the exact density of the system, including all non-linear ef- 
fects. To simulate laser induced phenomena it is customary to start from the 
ground-state of the system, which is then propagated under the influence of 
the potential 

VTD{r,t)=Ef{t)zsm{(jt). (114) 

vtd describes a laser of frequency uj and amplitude E 4 . The function f(t), 
typically a Gaussian or the square of a sinus, defines the temporal shape 



4 The amplitude is related to the laser intensity by the relation / = ^egcE 2 . 
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Fig. 4. Harmonic spectrum for He at A = 616 nm and I = 3.5 x 10 14 W/cm 2 . The 
squares represent experimental data taken from Ref. [65] normalized to the value 
of the 33rd harmonic of the calculated spectrum. Figure reproduced from Ref. [66] . 

of the laser pulse. From the time-dependent density it is then possible to 
calculate the photon spectrum using the relation 



where d(cu) is the Fourier transform of the time-dependent dipole of the sys- 
tem 



Other observables, such as the total ionization yield or the ATI spectrum, 
are much harder to calculate within TDDFT. Even though these observables 
(as all others) are functional of the density by virtue of the Runge-Gross 
theorem, the explicit functional dependence is unknown and has to be ap- 
proximated. 

5.2 High-harmonic generation 

If we shine a high- intensity laser onto an atom (or a molecule, or even a sur- 
face), an electron may absorb several photons and then return to its ground- 
state by emitting a single photon. The photon will have a frequency which 
is an integer multiple of the external laser frequency. This process, known as 
high-harmonic generation, has received a great deal of attention from both 
theorists and experimentalists. As the outgoing high-energy photons maintain 
a fairly high coherence, they can be used as a source for X-ray lasers. 

A typical high-harmonic spectrum is shown in Fig. 4 for the Helium atom. 
The squares represent experimental data taken from Ref. [65], and the solid 
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Fig. 5. Harmonic spectra of Hydrogen at a laser wavelength of A = 1064 nm for 
various laser intensities. Figure reproduced from Ref. [67]. 

line was obtained from a calculation using the EXX/EXX functional [66]. The 
spectrum consists of a series of peaks, first decreasing in amplitude and then 
reaching a plateau that extends to very high frequency. The peaks are placed 
at the odd multiples of the external laser frequency (the even multiples are 
dipolc forbidden by symmetry). We note that any approach based on pertur- 
bation theory would yield a harmonic spectrum that decays exponentially, i.e. 
such a theory could never reproduce the measured peak intensities. TDDFT, 
on the other hand, gives a quite satisfactory agreement with experiment. 

As mentioned before, high-harmonics can be used as a source of soft X- 
ray lasers. For such purpose, one tries to optimize the laser parameters, the 
frequency, intensity, etc., in order to increase the intensity of the emitted 
harmonics, and to extend the plateau the farthest possible. By performing 
"virtual experiments" , TDDFT can be once more used to tackle such impor- 
tant problem. As an illustration, we show, in Fig. 5, the result of irradiating 
a Hydrogen atom with lasers of the same frequency but with different in- 
tensities. For clarity, we only show the position of the peaks, and the points 
were connected by straight lines. As we increase the intensity of the laser, 
the amplitude of the harmonics also increases, until reaching a maximum 
at I = 1.5 x 10 14 W/cm 2 . A further increase of the intensity will, however, 
decrease the produced harmonics. This reflects the two competing processes 
that happen upon multiple absorption of photons: The electron can cither 
ionize, or fall back into the ground-state emitting a highly energetic photon. 
Beyond a certain threshold intensity the ionization channel begins to predom- 
inate, thereby reducing the production of harmonics. Other laser parameters, 
like the intensity, or the spectral composition of the laser, are also found to 
influence the generation of high-harmonics in atoms[66,67]. 
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5.3 Multi-photon ionization 




(a) (b) (c) 



Fig. 6. Ionization in strong laser fields: (a) Mult i- phot on ionization; (b) Tunneling; 
(c) Over the barrier. 



To better understand the process of ionization of an atom in strong laser 
fields, it is convenient to resort to a simple quasi-static picture. In Fig. 6 
we have depicted a one-electron atom at a time t after the beginning of the 
laser pulse. The dashed line represents the laser potential felt by the electron 
and the solid line the total (i.e. the nuclear plus the laser) potential. Three 
different regimes of ionizations are governed by the Keldish parameter, 

7=|- (H7) 

At low intensities (I < 10 14 W/cm 2 , 7>1) the electron has to absorb several 
photons before leaving the atom. This is the so-called multi-photon ionization 
regime. At higher intensities (/ < 10 15 W/cm 2 , 7 « 1) we enter the tunneling 
regime. If we further increase the strength of the laser field (I > 10 16 W/cm 2 , 
7< 1), then the electron can simply pass over the barrier. 

The measured energy spectrum of the outgoing photo-electrons is called 
the above threshold ionization (ATI) spectrum [68]. As the electron can absorb 
more photons than the necessary to escape the atom, an ATI spectrum will 
consist of a sequence of equally spaced peaks at energies 

E=(n + s)uj-I p (118) 

where n is a natural integer, s is the minimum integer such that sw-/ p > 0, 
and I p denotes the ionization potential of the system. 

Another interesting observable is the number of outgoing charged atoms 
as a function of the laser intensity. The two sets of points in Fig. 7 repre- 
sent the yield of singly ionized and doubly ionized Helium. The solid curve 
on the right is the result of a calculation assuming a sequential mechanism 
for the double ionization of Helium, i.e., the He 2+ is generated by first re- 
moving one electron from He, and then a second from He + . Strikingly, this 
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Fig. 7. Measured He + and He 2+ yields as a function of the laser intensity. The 
solid curve on the right is the calculated sequential He 2+ yield. Figure reproduced 
from Ref. [69]. 

naive sequential mechanism is wrong by six orders of magnitude for some 
intensities. 

Similar experimental results were found for a variety of molecules. Fur- 
thermore, in these more complex systems, the coupling of the nuclear and 
the electronic degrees of freedom gives rise to new physical phenomena. As 
an illustrative example of such phenomena, we refer the so-called ionization 
induced Coulomb explosion [70]. 

5.4 Ionization yields from TDDFT 

It is apparent from Fig. 7 that a simple sequential mechanism is insufficient 
to describe the double ionization of Helium. In this section we will show how 
one can try to go beyond this simple picture with the use of TDDFT [71]. 

To calculate the Helium yields we invoke a geometrical picture of ion- 
ization. We divide the three-dimensional space, It 3 , into a (large) box, A, 
containing the Helium atom, and its complement, B = 1R 3 \A Normaliza- 
tion of the (two-body) wave function of the Helium atom, \P(ri, r 2 ,t), then 
implies 




(119) 
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d 3 r 1 d 2 r 2 W{r 1 ,r 2 ,t)\ 2 , 
bJb 

where the subscript "X" has the meaning that the space integral is only over 
region X. A long time after the end of the laser excitation, we expect that 
all ionized electrons are in region B. This implies that the first term in the 
right-hand side of Eq. (119) measures the probability that an electron remains 
close to the nucleus; Similarly, the second term is equal to the probability of 
finding an electron in region A and simultaneously another electron far from 
the nucleus, in region B. This is interpreted as single ionization; Likewise, the 
final term is interpreted as the probability for double ionization. Accordingly, 
we will refer to these terms as p(°\t), p( +1 \t), and p^ +2 \t). 

To this point of the derivation we have utilized the many-body wave- 
function to define the ionization probabilities. Our goal is however to con- 
struct a density functional. For that purpose we introduce the pair-correlation 
function 

and rewrite 

P {0) (t) = \J Jd 3 r 1 d 3 r 2 n(r 1 ,t)n(r 2 ,t)g[n](r 1 ,r 2 ,t) 
p (+1 \t)= fd 3 rn(r,t)-2p(°\t) (121) 

J A 

p( +2 \t) = l-p(°)(i)-p( +1 )(i). 

We recall that by virtue of the Runge-Gross theorem g is a functional of the 
time-dependent density. Separating g into an exchange part (which is simply 
1/2 for a two electron system) and a correlation part, 

g[n]{r 1 ,r 2 ,t) = ^+g c [n]{r 1 ,r 2 ,t) (122) 
we can cast Eq. (121) into the form 

p<°>(t) - [N ls (t)} 2 + K(t) 

p {+1 \t) = 2JVi B (t) [1 - N ls (t)\ - 2K(t) (123) 
p (+2 Ht) = [1 - N ls (t)f + K(t) 

with the definitions 

N ls (t) = ^Jd 3 rn{r,t) = Jd 3 r |^ ls (r-,i)| 2 (124) 
K(t) =^JJd 3 r 1 d 3 r 2 n(r 1 ,t)n(r 2 ,t)g c [n}(r 1 ,r 2 ,t). (125) 
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Fig. 8. Calculated double-ionization probabilities from the ground-state of Helium 
irradiated by a 16 fs, 780 nm laser pulse for different choices of the time-dependent 
xc potentials. Figure reproduced from Ref. [71]. 

In Fig. 8 we depict the probability for double ionization of Helium calcu- 
lated from Eq. (123) by neglecting the correlation part of g. It is clear that 
all functional tested yield a significant improvement over the simple sequen- 
tial model. Due to the incorrect asymptotic behavior of the ALDA potential, 
the ALDA overestimates ionization: The outermost electron of Helium is not 
sufficiently bound and ionizes too easily. 

To compare the TDDFT results with experiment it is preferable to look 
at the ratio of double- to single-ionization yields. This simple procedure elim- 
inates the experimental error in determining the absolute yields. Clearly all 
TDDFT results presented in Fig. 9 are of very low quality, sometimes wrong 
by two orders of magnitude. We note that two approximations are involved 
in the calculation: The time-dependent xc potential used to propagate the 
Kohn-Sham equations, and the neglect of the correlation part of the pair- 
correlation function. By using a one-dimensional Helium model, Lappas and 
van Leeuwen were able to prove that even the simplest approximation for 
g was able to reproduce the knee structure [72]. As neither of the TDDFT 
calculations depicted in Fig. 9 show the knee structure, it seems that the 
approximation for the time-dependent xc potential is the most important in 
obtaining the ionization yields. 
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Fig. 9. Comparison of the ratios of double- to single-ionization probability calcu- 
lated for different choices of the time-dependent xc potential. Figure reproduced 
from Rcf. [71]. 

6 Conclusion 

In this chapter we tried to give a brief, yet pedagogical overview of TDDFT, 
from its mathematical foundations - the Runge-Gross theorem and the time- 
dependent Kohn-Sham scheme - to some of its applications, both in the linear 
and in the non-linear regimes. In the linear regime, TDDFT has become the 
standard tool to calculate excitation energies within DFT, and is now incor- 
porated into the major quantum-chemistry codes. In the non-linear regime, 
TDDFT is able to describe extremely non- linear effects, like high-harmonic 
generation, or multi-photon ionization. Unfortunately, some problems, like 
the knee structure in the yield of doubly ionized Helium, are still beyond the 
reach of modern time-dependent xc potentials. In our opinion, we should not 
dismiss these problems as failures of TDDFT, but as a challenge to the next 
generation of "density-functionalists" , in their quest for better approxima- 
tions to the elusive xc potential. 
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